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ABSTRACT 

We present a 3D numerical simulation of the recently discovered cometary structure produced as 
Mira travels through the galactic ISM. In our simulation, we consider that Mira ejects a steady, 
latitude-dependent wind, which interacts with a homogeneous, streaming environment. The axisym- 
metry of the problem is broken by the lack of alignment between the direction of the relative motion 
of the environment and the polar axis of the latitude-dependent wind. With this model, we are able to 
produce a cometary head with a "double bow shock" which agrees well with the structure of the head 
of Mira's comet. We therefore conclude that a time-dependence in the ejected wind is not required 
for reproducing the observed double bow shock. 

Subject headings: circumstellar matter - hydrodynamics ~ stars: AGB and post-AGB - stars: indi- 
vidual (Mira) - ISM: jets and outflows 



I 1. INTRODUCTION 

Martin et al. (2007) describe UV observations of a 
cometary structure extending over 2° in the sky, with its 
head centered on Mira. These observations were carried 
out with the GALEX satellite, with a filter with a cen- 
tral wavelength Ac = 1516 A and a FWHM AA = 256 A. 
Martin et al. (2007) speculate that part of the observed 
emission could be associated with the fluorescent UV 
lines of the II2 molecule. Also, part of the emission (par- 
ticularly in the head of the cometary structure) could 
correspond to emission from ionized species (e. g., from 
C IV) with lines which fall in the bandwidth of the fil- 
ter. We present a schematic diagram of the UV image of 
Mira's comet of Martin et al. (2007) in Figure 1. 

Mira is a binary system with an AGB star and a less 
luminous companion which is probably a white dwarf. 
The binary system moves through the plane of the galaxy 
with a spatial velocity of « 130 km s~^, at an orienta- 
tion of w 30° with respect to the plane of the sky, as 
can be deduced from the proper motion (Turon et al. 
1993), radial velocity (Evans 1967) and Hipparcos-based 
distance (107 pc, Knapp et al. 2003). The proper motion 
] is directed at an angle of 187° (measured E from N, see 
', Wareing et al. 2007). 

Wareing et al. (2007) computed a 3D simulation, in 
which they model this object as the interaction between 
an isotropic wind from the asymptotic giant branch 
(AGB) primary star and a streaming, homogeneous 
medium. Their simulation produces structures which are 
similar to the observed cometary structure. However, it 
is clear that the predicted flow does not show the de- 
tailed structure of condensations observed in the head of 
Mira's comet (see Figure 3 of Martin et al. 2007). 

In the present paper we explore whether or not the 
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"double bow shock" structure of the head of Mira's comet 
can be reproduced by a somewhat more complex model. 
In modelling the structures of planetary nebulae (PNe) , 
it is quite common to consider that the wind from AGB 
stars has a strong, equatorial density enhancement. Such 
a latitude-dependent wind can be used for reproducing 
the structures of many observed PNe and Proto-PNe 
(see, e.g., Balick & Frank 2002). 

In the case of Mira, we could also be seeing the interac- 
tion of a latitude-dependent wind with an environment 
which sweeps past the source. We have, however, few 
constraints of the characteristics and orientation of this 
latitude-dependent wind. The only real observational 
constraint is that the Mira binary system has an orbital 
plane that lies within ^ 30° from the line of sight, and at 
an angle of 130 — > 180° (measured E from N), with very 
high uncertainties (see Prieur et al. 2002). If we assume 
that the polar axis of the latitude-dependent wind lies 
perpendicular to the orbital plane, we would then have 
to place it close to the plane of the sky, and at an angle 
of 40 90° (measured E from N). This direction for the 
polar axis, taken together with the direction of Mira's 
proper motion (see above) implies that the polar axis of 
the outflow has to lie at an angle a ~ 33 — » 83° with 
respect to the direction of relative motion of the ISM 
which is streaming past Mira's wind. 

After exploring a limited set of ~ 10 parameter com- 
binations, we have chosen a model in which the angle 
between the polar axis of the latitude-dependent wind 
and the direction of the streaming environment has a 
value of 70°, which is consistent with the observations 
described above. In §2, we describe the choice of param- 
eters for our model, and the results from the numerical 
time-integration. In §3, we discuss the results. 

2. THE MODEL 

We model the cometary tail of Mira as the interaction 
of a latitude-dependent wind with a homogeneous envi- 
ronment in relative motion with respect to Mira. For the 
stellar wind, we consider a density of the form : 

p(r,0) = 4/(^), /(^?) = e-(e-i) |cos0|p, (1) 
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Fig. 1. — Schematic diagram of the UV image of the head of 
Mira's comet, showing the structures seen in Figure 3 of Martin et 
al. (2007). The proper motion of Mira is directed at an angle of 

7° to the W of S (N is up and E is left), approximately aligned 
with the large scale structure of the cometary tail. Therefore, when 
comparing this diagram with the model predictions (Figures 2-4), 
one has to allow for this small deviation from the vertical direction 
of Mira's motion. In this Figure, we introduce labels A-F in order 
to identify the different structures which are discussed in the text. 

where r is the spherical radius and 9 is the polar an- 
gle (measured out from the symmetry axis of the wind). 
In this equation, A is a scaling constant, ^ > 1 is the 
equator-to-pole density ratio, and p is a constant that 
determines the degree of flattening towards the equator 
of the density stratification. We have taken this form of 
the anisotropy function f{9) from Riera et al. (2005). 

We now consider a latitude-dependent wind velocity of 
the form 



v{e) 



(2) 



where vq is the velocity of the wind in the polar direction. 
With this choice, the ram pressure pv^ of the wind is 
isotropic. 

One can relate the scaling constant A with the mass 
loss rate of the wind AI^ through the relation 

p{r, e)v{r, 9) sm9d9 = 
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[f{9)]^^'^sin9d9. (3) 



We now choose an anisotropy function with p = 1/2. 
With this choice, one can compute analytically the inte- 
gral in equation ^ to obtain : 
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We now choose a ^ = 20 equator to pole density 
ratio, for which through equation [4] we obtain A — 
0.390M„/(47rwo). 

To complete the parameters for our anisotropic wind 
we choose a polar velocity vq = 10 km and a mass 
loss rate M„ = 7.7 x lO"'^M0yr"^ We set the wind tem- 
perature equal to lO'^ K at a distance ^ 8 x 10^^ cm 
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Fig. 2. — Volume rendition of the density structure obtained for 
a t = 8 X 10* yr time integration from the simulation described 
in the text. The horizontal extent of the displayed domain is of 



from the star, and impose the wind at all times within 
this radius. For the environment, we assume that it en- 
ters the grid along the z-axis, with a constant velocity 
= 130 km s~^, density rienv = 0.1 cm~'^ and temper- 
ature Tenv = 10^ K. 

Finally, we tilt the polar axis of the wind (out from 
which the angle 9 is measured, see equations [T] and [U 
so that it lies on an a;z-plane, pointing at an angle a — 
70° measured from the z-axis. This angle is consistent 
with the orientation of the orbital axis of the Mira binary 
system (see §1). 

We now carry out a 3D numerical simulation with the 
yguazii-a code in a (1.5, 1.5, 3.0) x 10^* cm domain. The 
stellar wind source is centered on an xy-plane, and placed 
at a zo = 6 X IQ-'^^ cm distance from the edge of the grid 
along the z-axis. This domain is resolved with a 5-level, 
binary adaptive grid giving 128 x 128 x 256 points at the 
maximum resolution. The yguazii-a code has been de- 
scribed in detail by Raga et al. (2000), and has been 
tested with laboratory experiments and employed for 
computing many different astrophysical flows (see, e.g., 
the review of Raga et al. 2006). We have used the version 
of the code in which the gasdynamic equations are inte- 
grated together with a rate equation for the ionization of 
H, and uses a parametrized cooling function (computed 
as a function of the density, temperature and H ioniza- 
tion fraction), which is described by Raga & Reipurth 
(2004). 

The numerical simulation is started with the wind (see 



Model for Mira's comet 



10" 



10' 



t=B,Oxia yr 



Bo 



10 5 10 
s [cm] 



-5. 10 5 10 
y [cm] 



10 5 10 
X [cm] 



Fig. 3. — Column density maps for t = 8 X 10* yr obtained by 
integrating the atom+ion column density along the y-axis (left) 
and along the x-axis (center). On the right we show the column 
density map for an orientation with the x-axis on the plane of the 
sky, and the 2-axis at an angle <p = 30° with respect to this plane. 
The column densities (in cm~-^) are shown with the logarithmic 
greyscale given by the bar on the top right. The origin of the 
coordinate system coincides with the (projected) position of the 
stellar wind source. 



equations [TJ2]) occupying all of the computational do- 
main, and the streaming environment entering from the 
bottom of the domain. Figure 2 shows a volume rendi- 
tion of the density of the flow after at — 8.0 x 10'' yr 
time-integration. The density stratification shows a lead- 
ing bow shock with a complex structure, and a tail with 
dense clumps and filaments. Figure 3 shows the column 
densities obtained from this stratification, calculated by 
integrating the density along the y- and the x-axes, and 
also assuming that the z-axis is at a = 30° angle with 
respect to the plane of the sky (with the x-axis on this 
plane). This orientation angle corresponds to the orien- 
tation of the motion of Mira with respect to the plane of 
the sky (see §1). The column density maps show struc- 
tures with large differences depending on the direction of 
the integration. 

In Figure 4, we show part of the time-evolution of the 
column density maps (calculated assuming a = 30° an- 
gle between the z-axis and the plane of the sky) of the 
region around the head of the cometary structure. In this 
time- sequence, we see that the region upwind from the 
star (i. e., below the star in the frames of Figure 4) at 
all times has a leading thick, dense region, which corre- 
sponds to the snowplow structure of the wind/streaming 
ISM interaction. The projection of the 3D structure on 
the plane of the sky has a complex shape, at some times 
giving a "double bow shock" morphology (see, e.g., the 
t = 7.0 X 10"* yr frame of Figure 4). This leading bow 
shock has a general side-to-side asymmetry that resem- 
bles the observations of Mira (see Figure 1). 

The separation between the star and the leading shock 
changes with time, giving values in the (3.6 4.2) x 
10^^ cm range (measured along the projected z-axis), 
corresponding to 3'. 8 — * 4'. 4 at the distance of Mira. An 
inspection of Figure 1 shows that this distance approxi- 
mately corresponds to the position of condensation A. 

The morphologies seen in some of the time frames of 
Figure 4 are remarkably similar to the head of Mira's 
comet. For example, the t = 6.5, 7.5 and 9.0 x 10* yr 
frames show a "double bow" head, in which one could 
associate one of the bow shaped structures with conden- 
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Fig. 4. — Time sequence of column density maps of the head of 
the cometary structure (the times in years are given on the bottom 
left of each of the frames). The column densities were computed 
for the </) = 30° orientation (between the z-axis and the plane of 
the sky, see the text and Figure 3) of the motion of Mira. The 
square domains which are shown have a size of 10^* cm, which 
approximately coincides with the horizontal extent of the diagram 
shown in Figure 1. The column densities (in cm~^) are shown with 
the logarithmic greyscale given by the bar on the top right. 

sation A, and the other one with condensation B (see 
Figure 1). 

In most of the time-frames shown in Figure 4, we see 
high density condensations in the region downstream (i. 
e., on top) of the star. These condensations could cor- 
respond to some of the knots labeled E in the schematic 
diagram of the observations of Mira's tail (Figure 1). 

The model which we are presenting produces a tail 
with a width of ~ 5 x 10^^ cm (in the (f> = 30° column 
density map, see Figure 3). This width is similar to the 
~ 7'. 5 (7.2 X 10^^ cm) width of Mira's cometary tail (see 
Figure 1 of Martin et al. 2007). 

3. CONCLUSIONS 

We have presented a 3D simulation of the interaction 
of a tilted, latitude-dependent wind with a streaming en- 
vironment. This simulation (viewed from an appropriate 
direction) produces column density maps which resemble 
the cometary structure around Mira in a qualitative way. 

Our simulation produces an asymmetric leading bow 
shock with a structure of condensations, and also a com- 
plex, time-dependent tail (see Figure 2). If one com- 
pares the predicted column densities with the UV maps 
of Martin et al. (2007), one can see that the individual 
time-frames (Figure 4) have condensations that fall at 
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the approximate positions of the observed condensations 
(see Figure 1 and Martin ot al. 2007). In particular, 
our model is successful at reproducing the "double bow 
shock" structure seen in the head of Mira's comet, with 
two bow shaped condensations at approximately the cor- 
rect distances away from the star. 

This comparison between the predictions and the ob- 
servations hinges on the the assumption that the emis- 
sion coefficient of the observed emission is approximately 
proportional to the density, so that the features seen in 
the column density maps can indeed be identified with 
the observed, emitting condensations. This assumption 
is of course likely to be incorrect. 

In order to advance in understanding Mira's comet it 
will be essential to obtain information about the mech- 
anisms responsible for the UV emission. Low resolu- 
tion spectroscopic observations will provide information 
about which mechanisms produce the emission of dif- 
ferent regions of the cometary structure. High resolution 
observations will provide radial velocity information, and 
proper motion measurements of the condensations of the 
cometary structure would also be possible by reimaging 
this object within the next few years. 

The parameters of our present model are not very well 
constrained, except for the known spatial velocity and 
the orientation of Mira's motion with respect to the sur- 
rounding medium. Also, if one assumes that the polar 
axis of the latitude-dependent wind is perpendicular to 
the orbital axis of the Mira binary system, one obtains 
an observational constraint for the orientation of the po- 
lar axis. We have chosen values for the mass loss rate 
and the wind velocity (and for the angular dependence 
of both of these) which are reasonable for an AGB star, 
and a value for the environmental density such as to give 
the right size for the head and tail of the cometary struc- 
ture. 

It is notable that through a very limited exploration 
(with ~ 10 attempts) of the choices for all of these pa- 
rameters, one is able to obtain a model that reproduces 
the remarkable "double bow shock" structure of Mira's 
comet. We therefore conclude that the observed dou- 
ble bow shock can be straightforwardly explained by a 
model of the interaction of a steady, latitude-dependent 
wind with a streaming cnviroimic;iit, and that it is not 
necessary to invoke the presence of a strong time depen- 
dence in the wind in order to reproduce the observed 
structures. Wareing et al. (2007) arrived at a similar 
conclusion from their analysis of the structure of the tail 
of Mira's comet. 

We end our discussion by mentioning the differences 
between our model and the model presented by Ware- 
ing et al. (2007). The main difference between the two 
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models is that while Wareing et al. (2007) considered an 
isotropic wind, we have considered the case of a latitude- 
dependent wind. Another, less important difference be- 
tween the two calculations is in the way the simulations 
were initialized. 

Wareing et al. (2007) started their simulation by im- 
posing a wind within a small, spherical region, having the 
rest of the computational domain filled in with a homo- 
geneous, streaming environment. This initialization can 
be seen as an idealization of having a star with a wind 
that suddenly "turns on" as it enters its AGB phase, 
but it is probably not appropriate for the case of Mira 
(which was probably in the AGB phase before starting 
to interact with the ISM in the galactic plane). 

We have started our simulation assuming that the 
stellar wind fills the computational domain, and that 
the streaming environment begins to enter the compu- 
tational domain from one of the grid boundaries. This 
is an idealization of the situation in which a stellar wind 
source that moves through a very low density environ- 
ment suddenly enters a much higher density medium (i. 
e., the star meets the galactic plane), which is also an 
idealized situation that probably does not correspond to 
the evolution of Mira's comet (which has been travelling 
through the stratified, galactic plane ISM for a much 
longer time than the time- integration of our simulation) . 

However, wc find that for integration times > 5 x 10* yr 
the structures obtained in the head of the cometary 
structure (see, e.g.. Figure 4) are qualitatively similar 
with both of the two initializations described above. A 
detailed description of the long tail of the cometary tail, 
however, might require a model in which the passage of 
Mira through a more realistic galactic plane ISM is con- 
sidered. 
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